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BOUNDED  ERROR  SCHEMES  FOR  THE  WAVE  EQUATION  ON  COMPLEX  DOMAINS 

SAUL  ABARBANELt,  ADI  DITKOWSKI*,  AND  AMIR  YEFET* 


Abstract.  This  paper  considers  the  application  of  the  method  of  boundary  penalty  terms  (“SAT”)  to 
the  numerical  solution  of  the  wave  equation  on  complex  shapes  with  Dirichlet  boundary  conditions.  A  theory 
is  developed,  in  a  semi-discrete  setting,  that  allows  the  use  of  a  Cartesian  grid  on  complex  geometries,  yet 
maintains  the  order  of  accuracy  with  only  a  linear  temporal  error-bound.  A  numerical  example,  involving 
the  solution  of  Maxwell’s  equations  inside  a  2-D  circular  wave-guide  demonstrates  the  efficacy  of  this  method 
in  comparison  to  others  (e.g.  the  staggered  Yee  scheme)  -  we  achieve  a  decrease  of  two  orders  of  magnitude 
in  the  level  of  the  Z^-error. 

Key  words.  Maxwell’s  equations,  wave  equation,  finite  difference-time  domain,  error  bounds,  boundary 
conditions,  complex  geometries,  staircasing 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Hyperbolic  systems  of  P.D.E.’s  describing  physical  situations  such  as  electro-magnetism, 
acoustics,  elastic  waves,  etc,  may  under  many  circumstances  be  cast  as  wave  equations  for  the  various  field 
components. 

One  class  of  problems  is  that  of  solving  numerically  the  Dirichlet  problem  on  complex  shapes,  e.g.,  inside 
wave  guides.  For  sufficiently  non-simple  geometries,  the  option  of  transforming  the  problem  to  body-fitted 
coordinates  is  not  always  a  viable  option,  especially  in  three  space  dimensions.  There  are  other  options,  such 
as  using  Cartesian  grids  and  approximating  the  body  shape  via  “staircasing” ,  “diagonal  split  cell  model” , 
etc  (see  for  example  Chapter  10  in  reference  [4]).  It  is  well  known  that  these  devices  are  not  very  efficacious, 
particularly  in  the  high  frequency  regime.  We  shall  demonstrate  that  “staircasing”  can  fail  even  for  low 
frequencies. 

In  this  paper  we  consider  the  application  of  the  method  of  boundary  penalty  terms  (“SAT” ,  see  references 
[1],  [2],  [3])  to  the  numerical  solution  of  the  wave  equation  in  a  finite  domain  with  Dirichlet  boundary 
conditions. 

In  Section  2  we  develop  the  theory  that  allows  us  to  use  a  Cartesian  grid  on  complex  geometries  and  yet 
maintain  the  order  accuracy  with  a  linear  temporal  error-bound. 

In  Section  3  we  construct  a  second  order  accurate  scheme  that  fulfills  the  conditions  imposed  by  the 
theory  presented  in  Section  2. 

Section  4  is  devoted  to  a  numerical  example  -  the  solution  of  the  transverse  magnetic  (TM)  Maxwell’s 
equations  [4]  between  two  concentric  circles.  (This  configuration  might  be  considered  as  a  cross-section  of  a 
very  long  wave-guide.)  This  problem  is  solved  using  four  different  numerical  algorithms.  Two  of  them  solve 
the  first  order  system  with  “staircasing”  -  the  Yee  staggered  scheme  [6]  and  a  4th  order  spatially  staggered 
scheme  due  to  Turkel  and  Yefet  [5].  The  other  two  solve  the  wave  equation  directly  on  a  non-staggered 
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Cartesian  grid,  one  with  the  SAT  formulation  and  one  without.  All  three  “standard”  (non-SAT)  algorithms 
have  very  large  errors;  the  SAT  algorithm  has  errors  that  are  at  least  two  order  of  magnitude  smaller. 
Summary  and  conclusions,  and  ideas  for  future  work  are  presented  in  Section  5. 

2.  Theoretical  Framework  of  the  Method.  In  reference  [1],  [2]  and  [3],  it  was  shown  how  the  case 
of  a  one-dimensional  P.D.E.  can  be  used  as  a  building  block  for  the  multidimensional  case  for  constructing 
error-bounded  algorithms  over  complex  geometries  with  Dirichlet  boundary  condition.  We  therefore  start 
with  the  following  one- dimensional  problem: 

(2-1)  S  =  S+/(M);  *>° 

(2.0a)  u(x,  0)  =  uq(x) 

d 

(2.0b)  —u(x,0)  =  uto(x) 

(2.0c)  u(rL,t)  =gi(t) 

(2.0d)  u(TR,t)  =  gR(t) 

and  f(x,t)  €  C2. 

Let  us  discretize  (2.1)  spatially  on  the  uniform  grid  presented  in  Figure  2.1.  Note  that  the  boundary 
points  do  not  necessarily  coincide  with  x\  and  xjy.  Set  Xj+ 1  —  Xj  =  /i,  1  <  j  <  N  —  1;  x\  —  Fl  =  7 xA 
0  <  7l  <  1;  -  xN  -  7 Rh,  0  <  7^  <  1. 


Fig.  2.1.  One  dimensional  grid. 


Since,  unlike  the  cases  discussed  in  [1],  [2],  equation  (2.1)  has  a  second  time  derivative,  attempts  to 
apply  naively  the  methods  presented  there  fail.  The  reason  is  that  if  we  follow  the  procedure  used  there  and 
write  the  following  discrete  approximation  to  (2.1), 

A 2 

(2.2)  ^-u  =  Du  +  f  (t)  +  Te 


where  u  is  the  projection  of  the  exact  solution  u(x,t)  onto  the  grid,  i.e.  u(xj ,  t)  =  Uj(t)  =  u(t);  and  write 
the  numerical  scheme 


(2.3) 


=  [Dv  -  tl{Alv  -  gL)  -  tr(Arv  -  gfl)]  +  f (t)  , 
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then  the  equation  for  the  error  vector  e  —  u  —  v  becomes 

(2.4)  7F-M‘  +  T- 
In  the  above,  v  is  the  numerical  approximation  to  u,  and 

(2.5)  M  —  D  —  tlAl  -  trAr  . 

D  is  a  differentiation  matrix  of  the  proper  order  of  accuracy  that  does  not  use  boundary  values.  The  matrices 
Al  and  Ar  are  defined  by  the  relations 

(2.6)  Alu  =  gL  -  Tl,  Aru  =  gR  -  Tr  , 

i.e.,  each  row  in  Al(Ar)  is  composed  of  the  coefficients  extrapolating  u  to  its  boundary  value  )  at 

TlO^r)  to  within  the  order  of  accuracy.  (The  error  is  then  Tl(Tr).)  The  diagonal  matrices  rL  and  tr  are 
given  by 

tl  —  diag(r//1 ,  7X2  j  *  *  *  *  )>  tr  —  diag (tr^  ,  tr2  ,  •  •  • ,  trn  )  . 

The  constrain  on  the  construction  of  the  A%  t’s  and  D  is  that  M  in  (2.4)  be  negative  definite.  The  negative 
definiteness  of  M  is  a  necessary  condition  for  extending  the  1-D  theory  to  the  multidimensional  case  (see 
[1],[3]).  Also  in  (2.4) 

(2.7)  T  =  Te  —  tlTl  -  trTr  =  (Ta,T2,  •  •  •  ,Tm,  ■  ■  ■  ,Tn)t  . 

If  the  matrix  M  can  be  diagonalized*,  then 

(2.8)  M  =  Q-1AQ 

with  the  diagonal  matrix,  A,  having  the  eigenvalues  of  M.  Defining  /x  =  Q  e,  equation  (2.4)  becomes 

d2  u  . 

^  =  Af.+<?r 

(2.9)  =  A  n  +  T  . 

This  is  an  un-coupled  system  of  O.D.E’s.  The  general  solution  for  the  mth  equation  is: 

limit)  =  cmieVI":t  +  +  -4=  [  sinh  (y/\ ~(t  -  s))Tm(s)ds  . 

V^m  JO 

Recalling  that  at  t  =  0,  e  =  e*  =  0  (i.e.  /x  =  fit  =  0  at  t  =  0),  the  solution  of  (2.9)  becomes: 

(2.10)  limit)  =  ~7f=  /  Tm(s)  sinh  [\/X^(t  -  s)]ds  . 

Note  that  unless  all  the  eigenvalues  of  M  are  real  and  non-positive  some  of  the  v^A^’s  will  have  a  positive 
real  part,  in  which  that  case  at  least  one  of  the  /zm’s  may  grow  exponentially  in  time.  In  order  to  prevent 
this,  we  have  to  demand  that  M,  in  addition  to  being  negative  definite,  also  possess  only  real  eigenvalues. 

*  Extensive  numerical  evidence  has  shown  that  theM  in  [1],[2]  (i.e.  representing  the  second  derivative  to^  and  2nd  order 
accuracy,  respectively)  has  distinct  eigenvalues  and  hence  is  diagonalizable. 
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Furthermore,  in  order  to  use  the  1-D  scheme  as  a  building  block  for  multidimensional  schemes,  M  should 
be  built  in  a  way  that  verifies  that  the  property  of  real  negative  eigenvalues  carries  over  to  the  multi¬ 
dimensional  differentiating  matrix.  One  way  to  achieve  this  goal  is  to  construct  M  as  a  negative-definite 
symmetric  matrix.  Then  an  estimate  on  the  error  bound  can  be  derived  directly  from  the  solution  (2.10), 


\l^m(t)\  < 


_ 1_ 

y/\\n 


LmM 


t 


where  TrriM  =  max0<s<(  |rm(s)|.  Then,  for  a  normalized  Q, 

(2.11)  ||  e||  =  ||  Mil  <  ^IITmP  , 


where  c0  =  minm=i,...,Ar  \/|Am|.  Therefore  ||  e||  grows  at  most  linearly  with  t. 

This  result,  of  a  linear  temporal  bound  on  the  error-norm,  can  also  be  derived  by  resorting  to  energy 
method  (see  [3]),  instead  of  directly  from  the  solution. 

Also,  as  mentioned  before,  the  construction  of  multi-dimensional  case 


d2u 

W 


—  V2u  +  /(x,t) 


on  complex  shapes  is  completely  analogous  to  the  method  indicated  in  [1],  [3]. 


3.  Construction  of  the  Scheme.  This  section  is  devoted  to  the  task  of  constructing  a  symmetric 
negative  definite  matrix  M  for  the  case  of  a  second  order  accurate  finite  difference  algorithm. 

Let 

'[1-210 
1-210 
0  1-21 
001-21 


1-2  1  00 
1-2  10 
1  -2  1 

1  -2  1 


'  0 

1 - 

O 

o 

o 

o 

1 _ 

C2 

1-3  3-1 

C3 

-1  4-6  4  -1 

cat-2 

-1  4-6  4  -1 

cat-i 

-1  3-3  1 

0  _ 

0  0  0  0  _ 
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(3.1) 


where 

(3.2) 
and 

(3.3) 


o 

o 

o 

> 

0  1-21 

-1  2  0-2  1 

-1  2  0-21 

-1  2-10 

0  0  0  0  _ 

j 

Cn-  1  -  C2  . 

Cfe  =  c 2  +  -  2)  , 


Cjy-I  -  C2 

N-3 


Note,  that  as  in  [2]  and  [3],  we  had  to  resort  to  using  connectivity  terms,  the  last  two  matrices  in  (3.1). 


(3.4) 


Al  = 


^(2  +  7l)(1 +  7l)  -7l(2  +  7l)  \{il  +  1l)  0  •••  0 


^(2  +  7l)(1  +  7l)  -7z.(2  +  7l)  0  •••  0 


(3.5) 


Ar  = 


o  ...  o  ^(7r  +  7r)  ~7r(2  +  7r)  ^(2  +  7r)(1  +  7r) 


0  ...  0  ^(7r  +  7r)  ~7r(2  +  7r)  ^(2  +  7b)(1  +  7r) 


(3.6) 


tl  =  -^2  diag  [rx,! ,  t£2  ,  tl3 , 0, . . . ,  0, 0] ; 


(3.7) 


TR  —  ^2  [^>  0,  . .  . ,  0,  2>  i )  TRw]  i 


In  order  to  make  the  matrix  M  =  D  —  tlAl  -  tr-A^  symmetric  we  choose: 


C2 


cat-i 


Tl2 


(1  -  7l)  7l 
2 

(1  ~  7r)  7r 
2 

3  -  Jl  ~  2  7l  rr,! 
1  +7i 
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(3.8) 


(3.9) 


T£/3 


tRn.  i 


TRn.  2 


-2  +  7 L±  1L  TLi 

2  +  t l 

TRn 

1  +  7 R 

-2  +  Jr  +  7 R  Tfijv 

2  +  7fi 


Tii »  >  4  • 


The  proof  that  the  symmetric  matrix  M  is  indeed  negative-definite  is  given  in  the  Appendix  to  this 
paper. 

Note  also  that  instead  of  solving  (2.3)  directly  as  a  2nd  order  O.D.E.  system  in  time,  one  can  solve 

^  =  [£>v  -  tl{Al\  -  gi)  -  tr(Arv  -  gfl)]  +  f 
dv 

—  =  w  . 
dt 

(3.10) 


The  number  of  ’variables’  has  increased  from  N  to  2N  but  one  gains  in  the  simplicity  of  the  time  integration. 


4.  Numerical  Example.  We  consider  the  dimensionless  Maxwell’s  equation  for  transverse  magnetic 
field  (TM,  see  [4])  in  two  space  dimensions: 


(4.1) 

(4.2) 

(4.3) 


d£_dHy_  dHx 
dt  dx  dy 
dHx  _  _dE 
dt  dy 
dHy  _  dE 
dt  dx 


where  Hx  and  Hy  are  the  x  and  y  components  of  the  magnetic  vector,  H,  and  E  is  the  electric  field  in  the 
^-direction.  The  set  (4.1)-(4.3)  is  to  be  solved  in  the  space  between  two  concentric  circles,  jjr  <  r  <  We 
consider  the  case  of  perfectly  conducting  boundaries.  Thus  the  boundary  conditions  are  given  by 


(4.4)  E&0,t)  =  0 

(4.5)  E(l,e,t)=  0. 

We  choose  the  following  initial  conditions  (note  the  polar  coordinates  r,  9): 


(4.6) 

(4.7) 


(4.8) 


E(r,  9,0)  =  cos#  [Ji(wr)  +  a  li(wr)] 

Hy(r,6, 0)  =  -sin2#|2~  [•Ma’r)  +oU(wr)] 


Hx(r,9, 0)  = 


[Jo(wr)  -  J2(wr)  +  a  Y0(u>r)  -  a  Y^ujr)]  j 
cos2  6 


cur 
sin2# 


\J\{u>r)  +  a  Yi(wr)] 


[Jo(wr)  -  J2(wr)  +  a  Yo(ur)  -  a  T2(wr)] 


where  the  Jn' s  and  the  Yn's  are  Bessel  functions  of  the  first  and  second  kind  of  order  n,  respectively.  Also, 


(4.9) 


a  Si  1.76368380110927;  u>  Si  9.813695999428405  . 
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The  exact  solution  of  the  IBV  problem  (4.1)-(4.8)  is  given  by: 

(4.10)  E(r,  0 ,  t)  =  cos (ut  +  0)  [J\{u)r)  +  a  Yi(oor)\ 

Hy(r ,  0,  t)  =  -  —  cos 0 cos(ut  +  0)  [Ji(ur)  +  a  Yi(wr)] 

(4.11)  +|  cos0sin(u;£  +  0)  [Jo(u>r)  -  J2(wr)  +  a  Yo(u>r)  -  a  ^(wr)] 

Hx(r ,  0, t)  =  —  cos  0  cos(cj£  +  0)  [Ji(o;r)  +  a  Yi(a;r)] 

(4.12)  — |  sin0sin(a;t  +  0)  [Jo  (or)  -  J2(wr)  +  a  Yo(ur)  -  a  Y2{u>r)\ 

We  note  that  we  can  extract  from  (4.1)-(4.3)  a  wave  equation  for  the  electric  field  E , 

d2E  _  a2£  d2E 
dt 2  0x2  (S^/2 

The  boundary  conditions  on  E  in  (4.13)  are  given  by  (4.4)-(4.5).  The  initial  condition  E(r,  0,  0)  is  given  by 
(4.6).  We  need  an  additional  initial  condition  on  Et)  which  we  obtain  by  differentiating  (4.10),  namely 

(4.14)  Et(r,  0, 0)  =  — a;sin0  [Ji(tor)  +  a  Yi(ur)]  . 


Four  numerical  schemes  were  used  to  solve  the  problem: 

(i)  The  Yee  scheme  [6].  This  second  order  accurate  scheme  is  staggered  both  in  space  and  time.  This 
entails  putting  initial  conditions  of  Hx  and  Hy  at  At/2  rather  than  at  t  =  0.  These  initial  conditions 
are  derived  from  the  exact  solution.  The  numerical  solution  is  carried  out  on  the  “staircased”  domain 
shown  in  Figure  4.1. 

(ii)  A  modification  of  the  Yee  scheme  (designated  Ty(2,4)),  see  [5].  This  one  has  4th  order  spatial 
accuracy  and  2nd  order  in  time.  The  stagger  and  the  “staircased”  domain  are  maintained  as  before. 

(iii)  The  SAT  algorithm  for  the  wave  equation  described  in  Sections  2  and  3.  The  grid  used  for  the 
numerical  integration  is  shown  in  the  right  side  of  Figure  4.1.  The  time  evolution  is  done  by  a  4th 
order  Runge-Kutta  method. 

(iv)  An  algorithm  which  formally  looks  like  the  SAT  in  (iii),  but  is  applied  to  the  “staircased”  domain 
of  Figure  4.1  (rather  than  SAT  one).  To  order  0(h2),  this  is  equivalent  to  using  a  standard  spatial 
central  differencing  scheme  with  the  nodal  points  at  edges  of  the  domain  given  the  boundary  value 
zero.  The  time  integration  is  done  as  in  the  case  (iii). 

We  first  present  the  L2  error  in  E  for  all  four  schemes  at  t  =  1  and  t  =  10  for  the  cases  Ax  =  Ay  —  h  = 
1/40,  h  =  1/80  and  h  =  1/160,  see  Table  1.  At  was  2/3  h  for  the  Yee  scheme,  h/ 18  for  the  Ty(2,4)  scheme 
and  hj 5  for  the  SAT  schemes. 

It  is  immediately  apparent  from  the  table  that  the  S AT-error  (scheme  iii)  is  at  least  2  orders  of  magnitude 
smaller  than  that  of  the  other  three  algorithms  at  all  the  various  times  and  grid  spacings. 

Since  the  non-SAT  schemes  have  errors  which  are  unacceptably  large  we  do  not  show  details  of  their 
temporal  behavior.  The  SAT  algorithm  (scheme  iii)  has  an  L2  error  which  grows  in  time  as  shown  in  Figure 
4.2.  We  see  that  this  temporal  growth  is  bound  by  a  linear  curve,  whose  slope  depends  on  h.  We  note  that 
for  all  reasonable  dimensionless  time  the  error  is  quite  small,  especially  for  h  <  1/80. 

5.  Conclusions  and  Discussion. 

(i)  It  seems  quite  clear  from  the  evidence  that  the  failure  of  the  non-SAT  schemes  is  due  to  the  fact 
that  “stair casing”  misrepresents  the  shape  of  the  body.  In  the  SAT  scheme,  on  the  other  hand,  the 
penalty  terms  take  account  of  the  true  shape. 
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Fig.  4.1.  The  “ staircased ”  domain  (left)  and  the  SAT  grid  (right),  h  =  1/40. 


A  =  1/40 

h  =  1/80 

h=  1/160 

t= 1 

i 

Yee 

0.4322 

0.3635 

0.1742 

ii 

Ty(2,4) 

0.4038 

0.3347 

0.1579 

iii 

SAT 

0.001203 

0.0001705 

1.5019e-05 

iv 

Staircased 

0.1022 

0.05041 

0.01936 

h=  1/40 

A  =  1/80 

A  =  1/160 

t=  10 

i 

Yee 

0.5101 

0.4364 

0.6683 

ii 

Ty(2,4) 

0.2642 

0.7079 

0.7243 

iii 

SAT 

0.008435 

0.0008354 

8.2707e-05 

iv 

Staircased 

0.7929 

0.4735 

0.7829 

Table  4.1 
The  Z/2  error. 


(ii)  The  numerical  results  validate  the  theoretical  predictions  of  the  temporal  behavior  of  the  L2  norm 
of  the  error. 

(iii)  Grosso-modo  the  CPU  time  per  node  is  of  the  same  order  for  all  schemes. 

(iv)  The  results  from  Table  1  and  Figure  4.2  seem  to  indicate  that  the  scheme  (iii)  converges  as  h 3, 
although  the  algorithm  has  a  truncation  error  of  order  h 2.  We  do  not  understand  this  pleasant 
anomaly,  although  it  is  possible  that  even  with  h  =  1/160  we  are  not  yet  in  the  asymptotic  conver¬ 
gence  regime. 

(v)  In  the  future,  we  would  like  to  apply  the  SAT  methodology  directly  to  hyperbolic  systems  such  as 
(4.1)-(4.3).  The  theory  is  not  complete  yet. 
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Fig.  4.2.  SAT,  L2  error  vs.  time. 


Appendix.  We  decompose  the  matrix  M,  defined  in  (2.5)  and  (3.1)  to  (3.8)  as  follows: 


) 


re: 


) 


) 


M  =  —  \olM\  +  (1  —  a)  M2  H~  M3  +  M4  +  M^] 

hr 


M2 


0  0  0 

0  0  0 

0  0-1  1 
1  -2  1 

1  -2  1 
0  1-100 
0  0  0 

0  0  0 


ml 1 

m4  2 

ml3 

ml2 

m2’2 

mf 

0 

(5.5) 

m4  = 

ml3 

2,3 

m4’ 

mf 

0 

0 

where: 


10 


0 


o 


(5.6) 

M5  = 

mN-2,N-2 

mN-l,N-2 

m^N~2 

0  mf-1’"-2 

mf-1’"-1 

«i r-1 

m^N~2 

N.N- 1 
m5 

rtf'* 

where: 


m 


£'N  =  l  +  2a- 


(1  +  7fi)  (2  +  7 r)  trn 


"5  -  *  .  —  2 
mN,N-i  =  -  a  +  jR(2  +  tk) trn 


mi 


n,n- 2  . 


7i?.  (1  +  Jr)  tRn 


m. 


2a  + 


77b  -  4  (l  +  7%,)  -  JR2  (2  +  77?)  (1  +  4 Tftjr) 
2 (1  +  tb) 


m 


N—l.N—2  _ 


=  1  —  a  — 


3t r  7ft 

2  +  2 


“I-  7ft  ^~ftjv 


N-2.N-2  o  ,  — 4  +  7b2  —  7R3  —  7fi2  (1  +  7ii)  tRn 

-2a+ - 2(2+73 - ' 

The  matrix  Mi  is  negative-definite  and  bounded  away  from  0  by  h2Tr2  by  the  argument  leading  to  eq. 
(2.4.31),  see  appendix  to  chapter  2  in  [3].  M2  is  non-positive  definite,  see  eq.  (2.4.34)  and  (2.4.35)  in  that 
appendix.  Prom  (3.2),  (3.3)  and  (3.8)  follows  that  ck  >  0,  k  =  1, . . .  ,N,  therefore,  the  matrix  M3  is  non¬ 
positive.  For  a  given  value  of  0  <  a  <  1,  tl,  and  trn  can  be  found  such  that  the  matrices  Mj  and  M5  will 
be  non-positive,  for  all  7 %  and  7b.  For  example:  for  a  =  1/10,  t^1  =  trn  =  4;  for  a  —  1/2,  tr1  =  trn  =  9 
and  for  a  -  8/10,  tLi  —  trn  =  24.  This  completes  the  proof  that  M  is  indeed  a  negative-definite  matrix, 
bounded  away  from  0  by  air2.  Therefore  the  norm  of  the  error  vector  ||  e  ||  can  grow  at  most  linearly  in 
time,  see  equation  (2.11). 
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